;+
; NAME:
;  hod_getmod
; PURPOSE:
;  Calls Takimitsu's fortran HOD model code with a given set of parameters
;  and returns a structure containing the output data
; CALLING SEQUENCE:
;  out=hod_getmod(modpar,modtab,[z=,cosmo1=,cosmo2=,rbin=,prof=,bmdl=,comm=,ver=])
; INPUT PARAMETERS:
;  modpar : n dimensional vector specificying N(M) model parameters
;  modtab : string specifying ascii table file name (in $HOD_DIR/tab/)
;           For HOD code v4.2a this is the parameter file.
; OPTIONAL KEYWORD PARAMETERS:
;       z : redshift
;  cosmo1 : vector of  h,OmeM,OmeL (default=[0.7,0.3,0.7])
;  cosmo2 : vector of sigma8_0,Ome_b,p (default=[0.85,0.045,0.97])
;    rbin : rmin,rmax,nbin (log binning/default=[0.05,80,60])
;    prof : 'nfw' (default) or 'm99'
;    bmdl : bias_model: SMT01, ST99 or TIN05 (default)
;    comm : command to call (default='wp_xi_ccf3_plth_tbl2')
;    ver  : hod code version (default=4.2a)
; OUTPUT PARAMETER:
;  structure with tags:
;    r_p,w_p_1halo,w_p_2halo,w_p_sum,xi_1halo,xi_2halo,xi_sum
; MODIFICATION HISTORY:
;  Created 1st Nov 2010 - J. Aird
;-
@hod_common
function hod_getmod,modpar,modtab,z=z,cosmo1=cosmo1,cosmo2=cosmo2,$
                    rbin=rbin,prof=prof,bmdl=bmdl,ver=ver,comm=comm

COMMON hod_model


if n_elements(ver) eq 0 then ver='v4.2a'
if n_elements(modtab) gt 0 then hod_parfile=modtab

;;;;;;;;;;;;;;;;;;;;;;;;
; V4.2a 
if ver eq 'v4.2a' then begin
   if n_elements(comm) eq 0 then comm='wp_xi_acf4_zheng07_i'
   ; all parameters are specified by the parameter file
   parfile=getenv("HOD_DIR")+"/parfiles/"+hod_parfile
   if not file_test(parfile) then hod_makeparfile,hod_parfile ; creates with default parameters
   
; calling hod code and reading output into a string array
   cd,getenv("HOD_DIR")+'/dat/HODdat' ; ideally don't want to have to move directories...
   command=getenv("HOD_DIR")+'/bin/'+comm+' '+$
           parfile+' '+$
           strjoin(strtrim(string(modpar),2),',')
   splog,command
   spawn,command,result,/stderr
   splog,'Completed command'
   print,result
   outstruct=hod_resultstruct(result)
   splog,'Read results'
endif
   

;;;;;;;;;;;;;;;;;;;;;;;;
; V3 code interface
if ver eq 'V3' then begin

   if n_elements(hod_modtab) gt 0 then begin
                                ; if the common block has been set use these values
      if n_elements(modtab) eq 0 then modtab=hod_modtab
      if n_elements(cosmo1) eq 0 then cosmo1=hod_cosmo1
      if n_elements(cosmo2) eq 0 then cosmo2=hod_cosmo2
      if n_elements(rbin) eq 0 then rbin=hod_rbin
      if n_elements(prof) eq 0 then prof=hod_prof
      if n_elements(bmdl) eq 0 then bmdl=hod_bmdl
      if n_elements(comm) eq 0 then comm=hod_comm
      if n_elements(ver) eq 0 then ver=hod_ver
   endif
   
; input checking and setting defaults
   if n_elements(z) eq 0 or n_elements(modpar) eq 0 or $
      n_elements(modtab) eq 0 then begin
      print,"USAGE: out=hod_getmod(z,modpar,modtab,[cosmo1=,cosmo2=,rbin=,prof=,bmdl=])"
      message,"Must specify input parameters z,modpar and modtab"
   endif

; checking specified modtab exists
   modtabfullpath=getenv("HOD_DIR")+"/tab/"+modtab
   if not file_test(modtabfullpath) then $
      message,"Cannot find file "+modtab
   
   if n_elements(cosmo1) eq 0 then cosmo1=[0.7,0.3,0.7] $
   else if n_elements(cosmo1) ne 3 then $
      message,"cosmo1 must have 3 elements"
   
   if n_elements(cosmo2) eq 0 then cosmo2=[0.85,0.045,0.97] $
   else if n_elements(cosmo2) ne 3 then $
      message,"cosmo2 must have 3 elements"
   
   if n_elements(rbin) eq 0 then rbin=[0.05,80,60] $
   else if n_elements(rbin) ne 3 then $
      message,"rbin must have 3 elements"
   
   if n_elements(prof) eq 0 then prof='nfw' $
   else if prof ne 'nfw' and prof ne 'm99' then $
      message,"prof must be 'nfw' or 'm99'"
   
   if n_elements(bmdl) eq 0 then bmdl='TIN05' $
   else if bmdl ne 'TIN05' and bmdl ne 'SMT01' and bmdl ne 'ST99' then $
      message,"bmdl must be 'TIN05', 'SMT01' or 'ST99'"
   
   if n_elements(comm) eq 0 then comm='wp_xi_ccf3_plth_tbl2'
   
; calling hod code and reading output into a string array
   cd,getenv("HOD_DIR")+'/tab'
   command=getenv("HOD_DIR")+'/bin/'+comm+$
           string(z)+' '+$
           strjoin(strtrim(string(modpar),2),',')+','+strtrim(modtab,2)+' '+$
           strjoin(strtrim(string(cosmo1),2),',')+' '+$
           strjoin(strtrim(string(cosmo2),2),',')+' '+$
           strjoin(strtrim([string(rbin[0]),string(rbin[1:2],format='(I5)')],2),',')+' '+$
           prof+' '+$
           bmdl
   print,command
   spawn,command,result
   
   outstruct=hod_resultstruct(result)

endif 






return,outstruct

end
  






    
 
